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_ ABSTRACT 

X 

H 

\ We have constructed realistic, self-consistent models of triaxial elliptical 

galaxies embedded in triaxial dark matter haloes. We examined three differ- 
ent models for the shape of the dark matter halo: (i) the same axis ratios as the 
luminous matter (0.7:0.86:1); (ii) a more prolate shape (0.5:0.66:1); (iii) a more 
oblate shape (0.7:0.93:1). The models were obtained by means of the standard 
orbital superposition technique introduced by Schwarzschild. Self-consistent so- 
lutions were found in each of the three cases. Chaotic orbits were found to be 
important in all of the models, and their presence was shown to imply a possible 
slow evolution of the shapes of the haloes. Our results demonstrate for the first 
time that triaxial dark matter haloes can co-exist with triaxial galaxies. 



Subject headings: galaxies: elliptical and lenticular, cD - galaxies: kinematics 
and dynamics - methods: numerical 
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1. Introduction 

According to the hierarchical picture of structure formation, galaxies originated from 
baryons that condensed within dark-matter haloes. If this picture is correct, dark matter 
should exist within and around galaxies. In the case of spiral galaxies, this has been widely 
confirmed by the study of rotation curves that remain approximately flat beyond the visible 
disk, suggesting the presence of invisible matter at large distances from the galactic centre 
(Rubin 1991). In the case of elliptical galaxies, the general absence of rotating gas disks 
makes this confirmation more difficult. Giant elliptical galaxies nevertheless show evidence 
of dark matter via stellar kinematics (Kronawitter et al. 2000), X-ray emission (Mathews 
& Brighenti 2003), and gravitational lensing (Keeton 2001). Numerical simulations of disk 
galaxy mergers including dark matter, like the recent work of Dekel et al. (2005), produce 
elliptical galaxies with stellar kinematics similar to what is observed. 

Large-scale simulations of structure formation typically exclude gas and stars, due to the 
computational complexities associated with gaseous dissipation, radiation, star formation, 
etc. These simulations make well-defined predictions about the structure of baryon-free 
dark matter haloes and the dependence of halo properties on mass and formation redshift 
(Navarro et al. 2004; Diemand, Moore & Stadel 2004). 

Much recent work has focussed on the shapes of the haloes (Allgood et al. 2006, Maccio 
et al. 2006; Bett et al. 2006). Simulated haloes are found to be generically prolate/triaxial. 
The mean short-to-long axis ratio is 0.5 ^ (c/a) ^ 0.7 and the distribution of c/a is ap- 
proximately Gaussian, with a sharp cutoff below c/a ~ 0.3 that probably reflects the onset 
of a dynamical bending instabifity (Merritt & Sellwood 1994). There are indications that 
mean halo elongation increases with halo mass (Allgood et al. 2006). The shapes of these 
simulated haloes are apparently supported mainly by anisotropic velocities rather than by 
flgure rotation (Warren et al. 1992). 

To the extent that luminous elliptical galaxies formed in a manner analogous to the 
(simulated) dark-matter haloes - i.e. through dissipationless clustering of smaller galaxies 

- one expects the structure of the two sorts of system to be similar. Gonfirmation of this 
was recently demonstrated in the case of radial density profiles (Merritt et al. 2005a,b): 
the same fitting function that best describes luminous elliptical galaxies is an equally good 
description of the A^-body haloes. Furthermore the simulated haloes appear to lie on the 
same relation between surface brightness and size (the "Kormendy relation" ) defined by the 
galaxies (Graham et al. 2006). 

With regard to shapes, observations of elliptical galaxies yield only the projected isophotes 
and so determination of the intrinsic shape distribution requires either statistical analyses 
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based on large samples (e.g. Ryden 1992, 1996; Tremblay & Merritt 1995, 1996; Ryden, 
Lauer & Postman 1993; Alam & Ryden 2002; Vincent & Ryden 2005; Plionis, Basilakos & 
Ragone-Figueroa 2006), or mapping of potentials via detailed kinematical data for individual 
galaxies (e.g. Franx, lllingworth & de Zeeuw 1991; Statler et al. 2004). Isophotal studies 
reveal that luminous and faint elliptical galaxies have significantly different shape distribu- 
tions: the latter have (c/a) ~ 0.75 while the former are flatter, (c/a) ~ 0.65 (Tremblay & 
Merritt 1996). The apparent shape distributions are found to be mildly inconsistent with 
axial symmetry in the sense that too few, nearly-round galaxies are observed, but these stud- 
ies do not strongly constrain the degree or distribution of triaxialities. However, detailed 
kinematical studies of individual galaxies sometimes find evidence for significant triaxiality 
(Statler et al. 2004). 

Still less is known about the true shapes of dark matter haloes, although constraints 
come from X-ray studies, lensing, and Milky Way kinematics (Tasitsiomi 2003). When 
comparing such data to the predictions of the iV-body models, it is important to keep in 
mind that halo shapes are likely to be influenced by the presence of the baryons. For 
instance, if baryonic matter accumulates at the center of the dark-matter potential, the halo 
is expected to become rounder and less triaxial (Merritt & Quinlan 1998; Springel, White & 
Hernquist 2004; Kazantzidis et al. 2004). 

Given these uncertainties, it is worthwhile to investigate new lines of evidence that 
might constrain the shapes of dark-matter haloes. Self-consistency studies constitute one 
such approach, which until now has not been explored. Can stationary, dark-matter haloes 
of arbitrary shape co-exist with a galaxy, or does the existence of a central, baryonic mass 
concentration imply restrictions on the elongation or triaxiality of the halo? 

In gravitational potentials that respect special symmetries (spherical- or axi-symmetry, 
for instance), isolating integrals of the motion exist in addition to the energy. Jeans's (1915) 
theorem provides a simple prescription for the construction of steady-state phase-space pop- 
ulations in such potentials: the phase-space density must be a function only of the isolating 
integrals. Orbits in triaxial potentials sometimes respect integrals in addition to the energy 
(Schwarzschild 1979), but the mathematical form of these integrals is typically unknown. 
Further, many orbits in triaxial potentials arc observed to respect no integrals aside from 
the energy (Merritt 1980; Valluri & Merritt 1998; Papaphilippou & Laskar 1998). 

Jeans's theorem can be generahzed to the non-integrable case: a steady-state galaxy 
is one in which the phase-space density is constant in every connected region (Kandrup 
1998). An example of a connected region is the invariant torus of a regular orbit, or the 
Arnold web corresponding to the chaotic parts of phase space at a given energy (Tabor 1989; 
Rasband 1990). Howevever Jeans's theorem says nothing about how this special state is to 
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be achieved, and stochastic orbits are often observed to remain confined to subsets of the 
full allowed region for long periods (e.g. Habib, Kandrup & Mahon 1997), making it unlikely 
that they would achieve a uniform phase-space population even in a Hubble time. 

In the current work, we use the method introduced in 1979 by Schwarzschild for the 
numerical construction of stationary galaxy models. Schwarzschild's method was tested on a 
modified Hubble profile (Schwarzschild 1979) and on a Plummer model (Siopis 1999). Later, 
it was used by Statler (1987) to demonstrate the self-consistency of model galaxies made 
of stars distributed according a triaxial density profile with a central core. Subsequently, 
high resolution observations of HST (Crane et al. 1993; Ferrarese et al. 1994; Lauer et 
al. 1995) showed that almost all the ellipticals have densities that rise toward the centre 
approximately as power laws p ~ r~"' , with 7 = 2 for fainter galaxies and 7 = 1 for 
brighter ones. Merritt & Fridman (1996; hereafter MF96) used Schwarzschild's method to 
construct self-consistent solutions galaxy models with both weak and strong cusps. In the 
case of strong cusps, the large number of stochastic orbits made it more difficult to obtain 
equilibrium models. Subsequent studies by Merritt (1997), Siopis (1999) and Terzic (2002) 
confirmed and extended the MF96 results. In these studies, attempts were made to represent 
the stochastic orbits via fully-mixed, steady-state distributions at each energy. 

All the works just cited dealt with galaxies containing a single (luminous) component 
only. Here, we introduce a second component, the dark matter, presumably consisting of 
long-lived, "cold," and coUisionless particles, e.g. supersymmetric neutralinos (Bertone, 
Hooper & Silk 2004). 

In the following Section we explain the method - a generalization of Schwarzschild's - for 
constructing self-consistent equilibrium models. The details of these models are discussed in 
Section [3] while in Section H] we describe their orbital characteristics. In Section O our results 
are analyzed and compared to those of models without dark matter. Finally in S ection E] 
we discuss the consequences of stochasticity for the long-term stability of the models. 

2. Method 

In the absence of analytic expressions for the orbital integrals, one must resort to nu- 
merical methods to construct self-consistent models of galaxies or dark-matter haloes. To 
obtain a numerical solution to the problem we follow the Schwarzschild (1979) scheme. The 
procedure can be summarized as follows: 

(1) A smooth, three-dimensional mass distribution resembling observed galaxies is pos- 
tulated and the corresponding gravitational potential calculated by means of Poisson's equa- 



5 



tion. 



(2) A library of orbits is realized by integration of the equations of motion in this 
potential starting from a large set of different initial conditions. 

(3) The configuration space is divided into cells and the time each orbit spends in any 
cell is recorded. 

(4) A linear combination of these orbits (with non-negative weights) is sought which 
gives the best discrete approximation to the known masses of the cells. 

We generalized Schwarzschild's method to the construction of two-component (luminous 
and dark matter), self-consistent models as follows. First, we constructed a catalog of orbits 
integrating the equations of motion for a particle moving in the potential generated by the 
sum of the luminous and dark matter density distributions. Then, two different grids of cells, 
one for the luminous and one for the dark matter components, were built. The technique 
used to realize each grid is the same described in MF96: the space was first divided by 21 
concentric ellipsoidal shells in zones containing an equal amount of one of the two matter 
components. Then, each octant was divided into three sectors corresponding to 

(i) the volume with ax > by and ax > cz, 

(ii) the volume with by > ax and by > cz, 

(iii) the volume with cz > ax and cz > by. 

Finally, we used a set of six planes to divide each sector in 16 regions. In this way, we 
obtained a total of 1008 cells per octant for each grid. The fraction of time spent by each 
orbit of this catalog in both sets of cells was recorded. We then used the data for the 
luminous-matter grid (called Luminous) to reproduce, by means of a linear superposition 
of orbits with non-negative occupation numbers, the given mass of the luminous component 
in each cell of that grid. The same procedure was followed to reproduce the dark matter 
distribution on the Dark matter grid. 

In practice, our optimization problem consisted of minimizing, separately, the quantities 




(1) 



and 
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where Bkj-ium is the fraction of time that the kth orbit spends in the jth cell of the Luminous 
grid and B^j-ctm is the same quantity in the jth cell of the Dark matter grid; Mj-ium is the 
mass which the model places in the jth cell of the Luminous grid and Mj-^^n that placed in 
the jth cell of the Dark matter one; Ck^ium and Ck-dm are the quantities to be determined by 
the minimization procedure. The latter represent the total mass of stars and dark matter 
particles, respectively, spread along the kth orbit (1 < ^ < norb)- 

To guarantee non-negative orbital weights, we imposed the constraints Ck-ium > and 
Ck;dm > 0. The NAG routine E04NCF (Stoer 1971; Gill et al. 1984) was used to solve these 
constrained least-squares problems. 



3. Mass models 

3.1. Density and potential 

We considered three different mass models for the galaxy+halo system. In each of the 
three cases, we adopted triaxial shapes for both components but with a different radial 
distribution for the dark and luminous matter. In the first model the dark matter was 
assumed to have the same axial ratios as the luminous matter, while the second and third 
models were characterized by haloes that were more prolate and more oblate, respectively, 
than the luminous matter. 

As in MF96, we adopted the following mass distribution for the luminous component: 

^ ^ M 1 

27iaibiCi m[l + my 

with 

x'^ xp' 

^'^ = — + + — 0<ci<bi<ai (4) 
"I 

and M the total luminous mass (see Fig. [1]) . 
For the dark matter component we adopted 

Pdm{m ) = — -— — (5) 

(1 + m'j (1 + m'^j 

where 

x"^ 

m''^ = ^— + To- + — < Cdm < bdm < adm, (6) 

^dm "dm ^dm 

and Pdmfl is the central dark matter density (see Fig. [T]). 
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The density profile in equation ([5]) was first proposed by Burkert (1995) for the dark 
matter haloes of dwarf galaxies, and later extended to the whole family of spiral (Salucci & 
Burkert 2000) and to massive elliptical (Borriello, Salucci & Danese 2003, hereafter BSD) 
galaxies. By adopting the Burkert profile - which has a large, low-density core - we are 
favoring data over theory, since A^-body simulations of gravitational clustering imply dark- 
matter densities that increase roughly as a power law inside of the halo virial radius (Navarro, 
Frenk & White 1996; Moore et al. 1998; Merritt et al. 2005b). Rotation curve studies, on 
the other hand, are generally interpreted as implying low (~ (1 — 5) x IO^^MqPc"^) dark 
matter densities within the region (r < 10^ pc) where the rotation curves are sampled (e.g. 
Burkert 1995; Salucci & Burkert 2000; de Blok & Bosma 2002; Gentile et al. 2005). It was to 
model these low observed dark matter densities that the Burkert model ([5]) was postulated. 

At large radii, the Burkert profile is similar to the NFW and Moore profiles that are 
commonly fit to A^-body haloes, i.e. pdm oc r~^. In the central regions where the two sorts of 
profile (core vs. power law) differ, the gravitational potential is dominated by the luminous 
galaxy. Thus, our conclusions about the constraints that a central galaxy places on the shape 
of the halo are probably applicable also to halo models without cores. 

The gravitational potential at a point x = (x, y, z) due to a mass distribution p = p{w), 
where w"^ = /a^ + y'^ /h"^ + /(? may be written (Chandrasekhar 1969) 





with 
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(9) 



Consequently, we find, for the luminous and dark matter potentials, respectively. 




(10) 



$dm(x) 





(11) 



with 




(12) 
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and 

A' 



' = \/(r + «L)(^ + &L)(^ + cL)- (13) 

The components of the gravitational forces are 



dxi Jo A(a^ + r)(l + m)3m 

OXi 

xf-^ (15) 

io A'(af + r)(l + m')(l+m'2) 

with oi = a;, 02 = 6/, 03 = q ,a'^ = a^^, «2 = ^dm; '^3 — ^dm and 2 = 1, 2, 3. These integrals 
were rewritten using an appropriate change of variables (see Appendix E]) in order to reduce 
the complexity of their numerical treatment. 

The mass contained in the ellipsoid of semi-axes wa, wb and wc is 

M{w) = A-nahc / p{w')w''^dw'] (16) 



^0 

this implies that the total mass is M^t = lim^^oo ^(w^)- In the case of luminous and dark 
matter, respectively, equation (fT6|l gives (see Fig. [1]): 

Mdm{m') = nadmbdmCdmPdm,oB{m'), (18) 

where B is the function 

B{m') = -2atan (m') + 2 log (1 + m') + log (l + m'^) . (19) 

We note that the total mass of the dark matter model is infinite. We cut off the dark matter 
grid at m = m'^^^, where m'^^g^^ ^19 was determined by the requirement that the mass in 
dark matter was ten times greater than the luminous total mass. 



3.2. Model parameters and constraints 

The free parameters of our models are: 



- ai, hi, ci and M for the luminous component; 
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- CLdm, bdm, Cdm and pdmfl for the dark matter. 



With regard to the axis ratios of the luminous component, there exist some constraints in 
the form of the observed distribution of galaxy isophotal shapes (Ryden 1992; Tremblay 
& Merritt 1995, 1996; Alam & Ryden 2002; Phonis, Basilakos & Ragone-Figueroa 2006). 
Intrinsic, short-to-long axis ratios of bright [Mb < —20) elliptical galaxies appear to be 
narrowly clustered around 3/4 (Tremblay & Merritt 1996); the constraints on the second axis 
ratio are less strong, but most studies (e.g. Tremblay & Merritt 1995) find that axisymmetric 
shapes can be securely ruled out. We accordingly fixed the short-to- long axis ratio of the 
luminous component to be 0.7 in all of the models and assumed "maximal triaxiality," i.e. 



^ = 0.7, T ^ = 1. (20) 

ai ai^ - ci^ 2 



For the dark matter, we investigated three choices for the axis ratios: equal to those 
of the luminous matter (MODI); a more prolate shape (T = 0.75) with Cdm/O'dm = 0.5 
(M0D2); and a more oblate shape (T = 0.25) with Cdm/o^dm = 0.7 (MODS) The parameters 
of the models are summarized in Table [TJ 

If the unit of length is taken to be a/, there remain two parameters: the ratio of luminous 
to dark-matter scale lengths and the normalization of the dark matter density. 

We applied the results of the fundamental plane study of Borriello, Salucci & Danese 
(2003) to the triaxial case, yielding adm/di = 3.64 and the value of pdmfl corresponding to 

^ MUK) ^ 3 (21) 
Mi (me) ^ ^ 

which is obtained by the evaluation of the luminous matter within the ellipsoid delimited by 
me, that characterizes the 'luminous' ellipsoid having the x-semiaxis equal to the effective 



TABLE 1 



MODEL 


Luminous matter 


Dark matter 




ci/ai 


Ti 


^dml ^dm '^dm 


MODI 


0.7 


0.5 


0.7 0.5 


M0D2 


0.7 


0.5 


0.5 0.75 


MODS 


0.7 


0.5 


0.7 0.25 



Table 1: Axial ratios and triaxial parameter of the density distributions of the 3 galactic 
models. 
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radius Re — 1.81a; (me is the triaxial generalization of the effective radius Rg), and that 
of the dark matter contained in the m'^ = 0.5 elhpsoid having the same long-axis extension 
as the luminous one. According to BSD this {adm/dh^e) paii' of values gives the best fit of 
theoretical data to the FP. From this constraint it follows that Pdmfl = 1-26 x 10'"^ Mi ^/ai^ 
for MODI, pdmfi = 2.26 x IQ-^Mi^T/ai^ for M0D2 and pdm,o = 1-14 x IQ-^Mi^r/ai^ for 
MODS. 

Hereafter we express all quantities in units such that G, ai and M are equal to unity. 
The implied unit of time is 

= G-'/W^'M-y^ = 1.49 X lOVrf^^r'^Y^V^' (22) 
' ^ VioiiMq/ Vlkpc/ ^ ' 



and Vu = \fGMjai is the unit of velocity; the, derived, energy unit is Eu = {GM'^/ai). 

4. Orbital catalogs 

4.1. Initial conditions 

The orbital catalogs were built as in Schwarzschild (1993). Initial conditions consisted 
either of zero initial velocity ("stationary" start-space), or points in the X — Z plane with 
Vx — Vz — and Vy — y^'2(E — $) (X — Z start-space). As argued by Schwarzschild 
(1993), these choices cover most of the orbital types in the full phase space of a nonrotating 
triaxial model; in particular, orbits starting in the X — Z plane are mostly tubes and avoid a 
region around the center of the model, while orbits starting on an equipotential surface are 
either stochastic or regular boxlcts, that approach the origin after a sufficiently long time. 
For each galactic model we fixed 32 energy values and for each value we selected a set of 
150 initial conditions from the X — Z start-space and a set of 192 initial conditions from the 
stationary start-space, thus obtaining a catalog of 10944 orbits. Each orbit was integrated 
in time using a 7/8 order Runge-Kutta-Fehlberg algorithm with a variable step size so as to 
have a relative error in energy less than 10~^ per time step. 



4.2. Integration times 

Particular care must be given to the choice of the orbital integration time. Because 
our aim is to construct a stationary model of a galaxy, all orbits must be integrated over 
a time interval long enough to guarantee a steady-state galaxy model. This means that 
the occupation numbers B^j of each orbit in the grid cells must be time- invariant. The 
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time required by the Bkj to converge to stationary values varies from one orbit to another, 
whether the orbit is regular or stochastic. In fact, a regular orbit can fill its invariant torus 
in a short time or, in the case of a nearly resonant orbit, it can take much longer to fill it 
(Merritt & Valluri 1999). In the same way, some chaotic orbits diffuse rapidly in the entire 
phase-space permitted by their energy, while others remain confined between regular tori for 
a very long time (Valluri &; Merritt 2000) (we will call ergodic the orbits with convergent 

To select the integration time for the orbits we followed an iterative procedure, proposed 
by Pfenniger (1984), that stops the integration when the relative variations of Bkj-jum and 
Bk,j-dm, in a time step go below 0.5%. In detail, we proceeded as follows: 

1) Integrate and store the Bkj-ium and B^j-dm over a time r = where is the 
"dynamical time," defined for each given energy as the period of the 1 : 1 resonant 
orbit in the X ~Y plane (see MF96); 

2) Restart the integration for another time r, and store the new B^j-ium and B^j-dm at 
the end; 

3) Compare the new and old Bkj-ium and Bkj-dm- If the maximum relative difference 
is greater than 0.005 then average the two sets of Bkj-ium and Bkj-dm, double r and 
repeat step 2, otherwise stop the computation. 

As argued by Pfenniger, this method does not converge for all orbits (e.g. "sticky" chaotic 
orbits), so it is needed to fix a maximum time of integration, tmax- 

We set this value at tmax = '^Th, where Th is the Hubble time; the choice Th = l5Gyr 
means tmax = 10,000Tu, where is the time unit in eq. [22] once assumed M = 10^^ 
and ai = Ikpc. 

Most of the CPU time was associated with the integration of the orbit libraries. Since 
each orbit is independent, this task is efficiently parallelized. We used the Plexus cluster at 
the Rochester Institute of Technology to carry out the integrations. The construction of the 
complete library of 10944 orbits for MODI required approximately 6 hours using in parallel 
16 nodes of the cluster. 
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5. 



Results 



5.1. Self-consistency 



The main result of this work is that we were able to succesfully construct self-consistent 
models in all of the cases considered for the shape of the DM halo. To obtain the model with 
the prolate shape of the dark matter distribution (MODS) we needed to enlarge the catalog 
of orbits to 11970. For each matter component, a parameter, 6, indicating the departure 
from self-consistency may be defined as 



where M represents the average mass contained in the grid cells and is the quantity 
defined in Section [2l The dependence of 6 on the choice of the number of orbits in the 
catalog is shown in Fig. [2] for M0D2; Fig. [3] gives the distribution of the relative error in 
mass, e, defined as e = \Mj — "^^=1 ^kBk,j\/Mj for each j'th grid cell in each of the three 
self-consistent models. 

If we analyze the orbits that contribute significantly to the self-consistent solutions, we 
find that a large fraction consist of orbits that do not conserve the sign of the components 
of the angular momentum (corresponding to regular box orbits and boxlets or to chaotic 
orbits). Regarding orbits that conserve the z— component (named short-axis tubes) or the 
X— component (outer tubes) of the angular momentum, it is noted that, to guarantee the 
prolate structure of the dark matter in M0D2, the fraction of outer orbits was greater than 
that of the short tubes, while the opposite was true in MODS, the more oblate model. 

To give a finer detail on the self-consistent model orbital structure. Fig. H] and [5] show 
the cumulative contribution, by mass, of the various orbital families in our self-consistent 
solutions. Chaotic orbits are of increasing importance at higher energies, especially for the 
dark matter component. This is partly explained by the initial abundance of high energy 
chaotic orbits in the full MODI catalog (see Fig. E]). Fig. E] indicates, also, how the low- 
energy back-bone of the self consistent model is mostly constituted by short-axis tube orbits; 
in the "more" evolved MODlbis model (described in more detail below), chaotic orbits of 
the luminous matter component are important on the whole energy range. All this suggests 
to investigate more deeply the role of orbital stochasticity, as we do in Sect. 6. 



S 




(23) 



M 



5.2. Velocity dispersion tensors 



To gain a more complete picture of the model kinematics, we also evaluated the first 
and second velocity moments in each of the self-consistent solutions. Using the properties of 
symmetry with respect to the coordinate planes in the cartesian frame, it was sufficient to 
store data for any given orbit in only one of the eight octants. 

We stored the velocities along an orbit symmetrizing them respect to the principal axes 
in order to insert the point in which one velocity is evaluated in a precise cell of the positive 
octant. In this way an orbit, starting with a fixed direction, correspond to a trajectory 
covered equally (or similar) in two opposite directions. It results a mean motions in each cell 
near to the zero value. Then we weighted these data with the occupation number (mass) of 
the luminous orbits in each cell obtaining: 



k=l ^ kdum^kj-^lum ^ '-'i ^ 



lum 



and 

k=l ^ k\lumJ-'k,j\lum - '-'i'-'h -^lum 
ViVh >lum= sr^Uorbit fi 75 • 

/ yfr— ] ^k;lurn-'-'k,j;lum 

Finally, we obtain the velocity dispersion tensor of the stars, crfum,ih ^^^h cell /, as: 

'^fumiih =< ViVh >lum,l ~ < Vi >lum,l< Vh >lum,l ■ 



The diagonalization of the velocity dispersion tensor, leads to the three principal velocity 
dispersions, crf^m,i;i — 1)2,3) having the profiles shown in the left column of Fig. [7] (we 
order the dispersions to have (jfum,i;i > '^Lm,/-2 > ^fum,i;3)- right column of the figure 

we show, for comparison, the corresponding principal dispersions evaluated in a model with 
only the luminous component. As expected, the velocity dispersion in the model with both 
luminous and dark matter have a fiatter profile. The difference in the external region is 
mainly due to insufficient orbital sampling and the (necessary) cut-off in the dark matter 
distribution. 



6. Importance of orbital stochasticity 

The Schwarzschild method for building self-consistent models uses orbits as "templates" , 
with supposedly time-invariant properties, so a requirement for a stationary solution is that 
a large fraction of the orbits (whether regular or chaotic) have nearly time-independent 
occupation numbers. Using the method of Pfenniger discussed above, we found that, at 
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the time 2Th-, roughly 20% of the occupation numbers of the orbits used to represent the 
luminous component were not yet fixed in time, while for the dark matter this percentage 
was as high as ~ 50%. For MODI, we attempted to obtain a more "stationary" solution by 
integrating the orbits up to 5Tf/. In this way, using a total of 2851 orbits for the luminous 
component and 2796 for the dark matter, we achieved a solution (MODlbis) with a relative 
error in the mass of each cell < 10~^ (see Fig. [8]). In this new model, the number of 
orbits with non-steady Bk^j values was reduced to ~ 15%, which would seem to indicate the 
presence in the initial solution of quite a few "sticky" chaotic orbits. 

Evaluating stochasticity based on orbital occupation numbers is prone to error, since 
even regular orbits can take long times to fill their invariant tori. To more accurately identify 
stochastic orbits, we used the Smaller Alignment Index (SALI) introduced by Skokos (2001), 
and later applied successfully to different dynamical systems (Skokos et al. 2003; Skokos et 
al. 2004; Manos & Athanassoula 2006). It is defined as (Skokos 2001) 



SALI{t) = min- 



Wi(t) 


W2(t) 


l|wi(t)|| 


I|w2(t)|| 



Wi(t) W2(t) 



Wil^tjll ||W2 

where wi(t) and W2(t) are two deviation vectors in phase space, centered in the initial 
condition of one orbit and pointing in two different arbitrary directions. In the case of 
Hamiltonian flows, SALI fluctuates around a non-zero value for regular orbits, while it tends 
(in time) to zero for chaotic orbits. 

By studying the evolution of this quantity for all the orbits of the catalog used to seek 
for the self-consistent solution of MODI we found that 68% of them are chaotic, especially at 
high energies (see Fig. [9]). Regarding the orbits of the subsample giving the self-consistent 
solution MODlbis, we found that more than 90% of non ergodic orbits are chaotic. 

To evaluate the degree to which the stochastic orbits in MODlbis would cause this model 
to evolve over long periods, we recorded the position and velocity of every chaotic orbit at 
the end {t = 5Th) of the integration and used these values as initial conditions for a new 



TABLE 2 



Set 


Luminous matter 


Dark matter 








X 


0.000872 0.00138 


0.00388 0.00137 


Y 


-0.00845 0.00150 


-0.01 0.00165 


Z 


0.00929 0.00177 


0.00771 0.00157 



Table 2: The mean values of the "evolution" parameter w and their standard deviations. 
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integration over another Th- We then again recorded the occupation numbers and compared 
them with the old ones. As shown in Fig. [TOl more than 90% of these "reintegrated" orbits 
have fluctuations in Bi j greater than the canonical value of 0.005. Operating in a similar way 
for dark matter, we reached the conclusion that our self-consistent model contains ~ 10% 
and ~ 23% of still-evolving stellar and dark matter orbits, respectively. 

In order to study the cell density change as a result of orbital evolution, and so the 
role of evolution on the overall shape of the model galaxy, we studied the variation of the 
occupation numbers of the chaotic-non ergodic orbits on a Th — 5Th time basis, of three set 
of cells (named X,Y,Z), sampling the volumes around the three coordinate axes. At this 
scope we defined, in each cell (i), an "evolution" parameter 



Wi = 







\lth 







(24) 



where in the numerator summation is extended over all the chaotic-non ergodic orbits in the 
self-consistent solutions and, to normalize, we divide for the sum over the whole set of orbits 
in the self-consistent solution at t = ^Th- The average value of w on the X, F, Z volumes gives 
an indication of the shape evolution: a globally unchanged shape would correspond to {w)x = 
{w)y = {w)z = 0. We found that for luminous matter < w >x is compatible with being = 0, 
while {w)y < and, {w)z > 0. This corresponds to a slight evolution toward a more prolate 
shape. The dark matter component shows a, statistically significant, increase in time of 
the chaotic orbit population of both the X and Z volumes, and a (smaller) depopulation of 
the Y volume. This should corresponds, too, toward a more prolate configuration. In view 
of future astrophysical applications (e.g. to the study of orbital the evolution of massive 
black holes or globular clusters subjected to dynamical friction), it is relevant to check how 
much the kinematical structure of the galactic models depend on the presence of a fraction 
of evolving orbits. To do this, we compared the principal velocity dispersions in MOD Ibis 
with those of MODI. As we can see by comparing Fig. [11] with Fig. [TJ the behaviours are 
similar, with a relative difference that is almost always below 10%. We checked that the 
largest relative differences refer to the outermost shell. 
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A. Potential and forces 

Using tlie substitution s = (1 + r)"^/^ in Eq. [TO]and[TT| they transforms into: 

•I'/Cx) = -GM 



ds 


1 


L ^ 


(1 + m)2_ 



(Al) 



where 



and 



^dmi^) = -'^TrGadmbdmCd 



'mPdmfl 



X J ^ — atanm' + log(l+m') — ^log(l+m'^^ 



m^(s)=s^ 



1 + S2(a;2 - 1) 1 + S^{bi^ - 1) ' 1 + s2(q2 - 1) 



A=J[1 + s\ai^ - 1)][1 + s^ik^ - 1)][1 + s2(q2 - 1)] 



X 



A'= J[l + S\ad^^ - 1)11 + _ ^ ^2(c,^2 _ 1| 



(A2) 



(A3) 



(A4) 



With the substitutions s = ai-i{al- + T) i^j^g luminous forces and s = adm;iio.dm;i + 

r)^^/2 for that of dark matter we transform Eq. [T4l and [TSl in: 



(9$, 

F;, = -GM-^ = -2x, 

OXi 



s^ds 



/^imiil + rrii 



(A5) 



drrii 



dx 



X 



4:7vG(ldmbdmCdmpdm,0-^i 

s^ds 

7^ 



AKl + m9(l + mn 



(A6) 
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where xi = X2 = y, xs = z and: 



mi^{s) = 



[s) = s 



^ 4. f j_ ^ 



Ti 



;3S 



The costants in Eq. IA7I are 



2 = 1: 

Aj-i = 6/ — 1 Aj.2 = c/ — 1 = 
~ ^ ^i;2 ~ ~ 1 C*j;3 = Cj^ — 1 

i = 2: 

^r,i = - Aj,2 = 1 - = 6/ 



i = 3 : 

Cj-]^ = 1 — Cj'^ ^j,2 ~ bj — Cj'^ ~ 0- 

with j = I 01 j = dm if we deal with luminous or dark matter, respectively. 



(A7) 



(A8) 
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Fig. 1. — Top: long- axis density profiles of the luminous and dark matter components. 
Middle: mass enclosed within spheres of radius r in the case a = h = c. Bottom: circular 
velocity profiles in the spherical model, including both components (dashed line), and in a 

mndfil withniit dark mattfir 
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Fig. 2.- 



The rms deviation 5 (equation [23|) as a function of the number of orbits supplied to 
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Fig. 3. — Histograms of the relative errors (see text). They refer to both the components 
grids of the three models analyzed and are calculated using the quantities of the self- 
consistent solutions. 
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Fig. 4. — Cumulative, by mass, energy distributions of the various orbital families in the 
self-consistent solution of MODI (left column) and MODlbis (right column). The symbols 
"B," "X," "Z," and "C" denote the mass contributed by box, x-tube, z-tube, and chaotic 
orbits, respectively. 
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Fig. 5. — Cumulative, by mass, energy distributions of the various orbital families in the 
self-consistent solution of M0D2 (left column) and MODS (right column). The symbols "B," 
"X," "Z," and "C" denote the mass contributed by box, x-tube, z-tube, and chaotic orbits, 
respectively. 
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Fig. 6. — Cumulative, by number, energy distribution of the whole MODI orbital catalog. 
Symbols as in Fig. H] 
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Fig. 7. — The stellar principal velocity dispersions crfumi i^ ^ — 1,2,3 as functions of the 
distance from the centre in model MODI (left column) and in a model with the same 
luminous mass but without dark matter (right column). 
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Fig. 8. — Histograms of the relative errors for MODlbis, referring to the luminous and 
dark matter component, respectively. 




Fig. 9. — Fraction of different types of orbits in the full MODI orbital catalog vs energy. 
Full circles represent stochastic orbits; open circles are box-like orbits; open and full squares 
correspond to short and outer tube orbits respectively. 
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Fig. 10. — Histograms of the maximum |ASj_fe| for the stochastic orbits used in MODlbis 
integrated once up to ST^ and once up to 67^. 
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Fig. 11. — Left column: principal velocity dispersions for the luminous component, obtained 
with the self-consistent solution realized after having integrated orbits over ^Th (MOD Ibis). 
Right column: relative change in this quantity with respect to that obtained with the self- 
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